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1. Introduction 

Parameter-dependent evolution of spectra in random matrix ensembles has been a 
subject of investigation almost from the inception of the Random Matrix Theory 
(RMT) Dyson pp showed that the joint probability distribution function of 

the eigenvalues in the classic random matrix ensembles can be alternatively obtained 
as an equilibrium distribution of a one-dimensional gas of classical charges which 
undergo Browninan motion in fictitious time while simultaneously logarithmically 
repelling each other. The same dynamics of eigenvalues arises if the corresponding 
random matrices themselves execute Browninan motion [3| according to a sequence 
of mappings H i— ► H + dV where H is drawn form a random matrix ensemble and 
dV is the elementary step of the random walk in the space of matrices of the same 
symmetry as H. 

However, despite the theoretical interest in these types of models, the 
experimental motivation in the early days of RMT was lacking since the nuclear 
spectra for which RMT was originally devised 0] are essentially immutable. Later, 
the range of RMT applications has broadened to include systems such as mesoscopic 
quantum structures (e.g., quantum dots) which can be subject to a variety of external 
perturbations In fact, controlled external perturbations such as magnetic field 
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or gate potential often serve as important experimental tools to study statistical 
properties of such systems. Thus, instead of a sequence of random steps dV one 
is lead to consider finite mappings H i— > W = H + V , where V is a fixed matrix with 
the same symmetry as the ensemble from which the matrices H are drawn f jS] [7| |S] . 

While magnetic field clearly belongs to the class of extended perturbations, in 
many applications the assumption that the rank of V scales with TV may not prove 
to be adequate in the RMT modelling of mesoscopic devices. External perturbations 
such as an STM tip, a change in the potential on a finger-shaped gate, or a jump in 
the position of a bi- or multi-stable defect represent potentials which are localised in 
space. In the RMT language localised potentials are modelled by matrices V of finite 
rank r. Such matrices can be represented using a finite set of orthonormal complex 
vectors {a K }^ =1 as V = Nj^^i v k.&-k ® a«- [For technical reasons related to the use 
of the Itzykson - Zuber - Harish-Chandra (IZHC) integral 10 , our consideration here 
is restricted to the unitary ensemble, hence the requirement that a K are complex.] 

In the standard RMT one can distinguish two classes of correlation functions 
between the levels 0- In the first class the levels are "labelled" by their positions 
("energies" in the Hamiltonian interpretation of the random matrices), so that the 
fundamental object is the many-point correlation function of the density of states 
(DoS)u(E)=ZLiHE-e a ), 

R n ({E}U) = (flHEiY), (1) 

where the angular brackets denote averaging over the matrix ensemble. In the bulk 
region of the spectrum the correlation functions are translationally invariant with 
respect to the simultaneous shifts of all energies, so that R n is actually a function of 
n — 1 energy differences. Such correlation functions are known to have determinantal 
structures in all three classic Wigner-Dyson ensembles. E.g., in the unitary ensemble, 
the N — > oo limit results in 

R n ({E}2 =1 )=detk{e i -e j ) (2) 
y 

where k (e, — Sj) = k (si — Sj) — 5 {en — £j) 9 (i — j), k (u) = suittu/ttu is the celebrated 
sine kernel of the Wigner-Dyson RMT, the dimensionless energies e = E/A are 
normalised by the mean level spacing A = (y (E)) 1 = p (E) , and the step-function 
9 is regularised so that 9(0) — 0. The term involving the <5-function accounts for 
the self-correlation of levels. Analogous structures with more involved kernels arise 
in other ensembles. While originally derived in the context of the Gaussian Unitary 
Ensemble (GUE), the correlation kernel k acquires the same universal form in the 
bulk scaling limit for a wide class of non-Gaussian random matrix ensembles 

V(H)=C lJ e- TlW W, (3) 

where W(H) is an even polynomial, and C\j is the normalisation constant. 

In the second class of correlation functions the levels are labelled by their numbers, 
while their energies may or may not enter as additional variables depending on the 
specific correlation function. An archetypal example is the distribution of level 
spacings &i(u>) which is a particular case of the distribution & q (to) of distances 
between a level a and another level separated from a by q — 1 other levels: 

& q (u)=A(8(e a -E)S(e a+q -E-u))). (4) 

f A separate subject which is beyond the scope of this paper concerns mappings such that the 
symmetry of V is lower than the symmetry of H . 
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©i is actually a correlation function between an arbitrary level a and its neighbour 
a+1. Translational invariance in the bulk of the spectrum ensures that this correlation 
function does not depend on either E or a. In contrast with no summation over 
a is involved, so that level-number-specific information is retained, and, indeed, near 
the edge of the spectrum, the same correlation function - now as a function of both 
level numbers - describes the distribution of distances between, say, the second lowest 
and the third lowest levels. 

The two types of the correlation functions are related to each other via infinite 
series. Thus, an equivalent expression for ©i(w) without referring to specific level 
numbers is the joint probability to find two levels at positions E and E + u> and all 
other levels in the interval M./[E,E + u]. After straightforward algebra 0, ©i(w) is 
written as 

SiM = -Kdl T / dEx... dE n R n {{E t } n i=l ) , (5) 

where i? n ({E}™ =1 ) = dety k (Ei — Ej) is the part of R n from which the level self- 
correlation terms are excluded. Combined with the determinantal structure of R n , 
this formula allows one to express ©i in terms of the derivative of a certain generating 
function, ©i = — det[B,E+u] (l — kj, where the determinant is understood as 
a Fredholm determinant on the space of functions with the support in the interval 
[E, E + u>], and the kernel k is defined through its matrix elements: k(u, w) = k(u—w). 
Analogous expressions can be obtained for other P q 's. The Fredholm determinant 
itself - the generating function - is the probability to find no levels in the interval 
[E, E + u>]. Other © g 's are similarly related to the Fredholm determinant expressions 
for the probability to find q—1 levels in a given interval. 

Conversely, the correlation function R2{E, E + oj) can be expressed as an infinite 
sum over © g 's, 

oo 

R 2 {E,E + u)= 6 ?H ( 6 ) 

q— — oo 

which is a simple restatement of the fact that the conditional (on the existence of a 
level a at E) probability to find any level at E+lu is the sum over q of the (conditional) 
probabilities to find the (a + q)th level at that position. Operationally, however, R n 
is usually more readily accessible than 6 q , so that equation 10 and its analogues are 
not often used in practical calculations. 

In direct analogy with the duality between the two classes of correlation functions 
described above for ordinary RMT, one can introduce their generalisations in the 
parametric RMT. Thus, instead of R n we are now lead to define the parametric many- 
point correlation function as 

In m \ 

R nm (mu, {e' 3 }?=i) = ( n v (^) n v> m) ) > ( ? ) 

\i=l j=l I 

where v' (E) = TtS (E - H'). 

Except in the r = 1 case (see below), the "level spacing" in the combined set 
{E} U {E 1 } does not have an immediate physical interpretation. A more meaningful 
analogue of equation is the distribution of distances between a level a and the 
parametric "descendant" of the level a + q, 

Pq (w, f) = A (6 (e a - E) 5 (e Q+9 (r) - E - u)) , (8) 
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where y is the set of the parameters characterising V. Due to universality 0^], y 
can be identified with either a single parameter x measuring the overall strength of 
V in the case of extended (infinite rank) perturbations, or with a set of at most r 
parameters in the finite rank case (see Section 3). Pq is simply the distribution of 
shifts of a single level under the influence of a parametric perturbation. Note the 
absence of summation over a; nevertheless, in the bulk of the spectrum P q is not a 
function of a. 

A related quantity which in the past has been extensively studied in the context 
of parametric correlations in chaotic and disordered systems using either numerical 
or approximate analytical techniques [131 1141 1151 116j is the single level velocity 
correlation function Co (x). Within the formalism developed in the present paper, it 
is convenient to represent it as a special case of the following correlation function: 

C 9 (r) = (d ? e Q ( y = 0)d f e Q+9 (r)), (9) 

which represents the correlation of responses to the perturbation of two levels 
separated by a distance y in the parameter space, and by q — 1 levels in the level 
sequence. 

The objective of this paper is to present a set of interrelated developments in the 
parametric RMT. In the next section we construct a formalism to determine the joint 
distribution functions for the eigenvalues of H and H' for perturbing matrices V of 
arbitrary rank r. In particular, we concentrate on the case when r is finite in the 
N — > oo limit, although the results are equally applicable in the opposite limit where 
we reproduce the Gaussian transition kernel implicit in the earlier work on the subject 

ei he]. 

In Section 3.1 we use the joint distribution functions to derive the many-point 
correlation function R nm . We show that the latter possesses the determinantal 
structure at arbitrary r |17| . generalizing the result previously known only in the 
context of two-matrix models |3 I18j (i.e., with V itself drawn from a random 
matrix ensemble), and thus corresponding to r — N in the classification adopted 
here. In Section 3.2 these results are applied to prove (on a physical level of 
rigour) that the universality of the Wigner-Dyson sine kernel extends to a class 
of non-invariant ensembles characterised by the distribution function Vy (H) = 
Cu exp{-Tr W(H -V)}. 

The final objective of this paper is to present (Section 3.3) a general framework 
for the calculation of P g (w,y), C q (y), and related number- dependent correlation 
functions in terms of the corresponding family of generating functions. Owing to 
the determinantal structure of R nm |17l ^] (see also Section 3.1 below), the 
latter possess representations in terms of Fredholm determinants of certain parameter- 
dependent integral kernels. 

2. Joint Distributions of Eigenvalues 

Consider a set of pairs of matrices H and H' = H + V, with H being a random 
Hermitian N x N matrix H drawn from the unitary ensemble V (H) and V is a 
deterministic matrix of rank r. The joint distribution of the matrices H and H' is 

V (H, H') — V (H) 6 (H + V - H') , (10) 

where the matrix delta-function is understood as the product of ordinary scalar delta- 
functions - one per each of the N 2 independent variables of H . The joint distribution 
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function of the combined set of eigenvalues {ea} a= i of the matrix H and j of 

H' is obtained from the above expression by integrating out the angular degrees of 
freedom: 

V ({6} , W}) = I A 2 (H) A 2 (H') dfi (U) dy (W) V (H, H') , (11) 

where dfi (U) is the invariant Haar measure on the unitary group U(N), and A (H) is 
the Vandermonde determinant 

A (H) = A ({e}) = [] i^-ep). (12) 

l<a<0<N 

With the help of a Lagrange multiplier matrix A (it is straightforward to show that 
independent components of A can be arranged into a Hermitian matrix) the expression 
for the joint distribution function takes the form 

V ({e} , {e'» = C V f A 2 (H) A 2 (H>) A 2 (A) dp (U) dy (V) ^ ^ f[ 

J (21") k=i 77 

x exp {i Tr (AH + AF - AH')} ex P {~W (H)} , (13) 

where A^ are the eigenvalues of A. Integration over the angular degrees of freedom of 
H and H' can be performed with the help of the IZHC integral [T(]| 

I a (tt\ l-rr attuttU det a/ 3 [exp (ia a bp)] 
J dfj, (U) exp {iTtAUBW } = c — — — — , (14) 

where {o a } a=1 , {bp}p =1 are the complete sets of the eigenvalues of the matrices A and 

B, respectively, and c = (i/2)( N _JV V 2 YljLi 3- The remaining integral thus takes the 
form 

V ({e} , {e'» = / A (H) A (H') d M (U x ) {[ d\ k exp {iTr (XV)} 

( 27r ) J fe=i 

x (-1) CT+CT exp ^ i ^ (e Q A (T ( Q ) - e' Q A^ (Q) ) I exp ^ - ^ W (e Q ) I , 

er,a-' L a=l ) \ a—1 ) 

where the sum over a runs over all permutations of the indices 1, ...,N, and ( — 1) CT 
denotes the signature of the permutation. The ordering of the energy levels in the 
Vandermonde determinants is assumed to be in the increasing order. Due to the 
invariance of the measure /j, (U), the integral in l|15fl does not depend on the "angular" 
degrees of freedom of V. The latter can thus be thought of as "gauge" degrees of 
freedom of l|15fl . and integrated over with an appropriate measure d/i( r ) ({a K }). A 
convenient choice for dfir r -\ ({a K }) is the uniform distribution of all components of 
a K constrained only by the requirement of orthonormality. U\ rotations can now 
be absorbed into the rotations of V, and the integral over dn(U\) gives a constant 

t L N = J dpb(JJ\) — (271")^ ~ JV )^ 2 /YijLiJ" Equation (|15fl can now be rewritten as 

V ({e} , {e'» = ^p^e-E2- (H) A (H') 

(2tt) 

x J d/X( r ) ({a K }) det S (e a - e' - V aa ) , (16) 



(15) 
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where V aa are the diagonal elements of V = NY^ K= i ® a t- ^ is worth noting 
at this point that if integrating over dfj,/ r ) ({a K }) does not introduce correlations 
between V aa , the distribution function acquires, apart from the asymmetry introduced 
by the confining potential, a form reminiscent of a quantum-mechanical transition 
amplitude: the Vandermonde determinants play the role of Slater determinants 
describing antisymmetric many-fermion wave functions, and the determinantal form of 
the transition kernel implies that the fermions are non-interacting albeit subject to an 
unusual single-particle Hamiltonian. Correlations between V aa then imply interactions 
between the particles. 

We now analyse several limiting^ cases. If V itself is drawn from a unitary random 
matrix ensemble V(V) oc e~ TlV l 2X (and is thus of rank r — N), the integration over 
cfyi(Af) ({a K }) is subsumed into the integration over V. In "Cartesian" coordinates 
Ha dV aa Y[ a <f3 dMVafid^Vap the matrix elements V aa are explicitly uncorrelated 
and independently Gaussian distributed. The transition kernel is thus proportional 

— X 2 (e -c' "| 2 /2 

to deta/3 e \ a ®> i corresponding to free propagation of the ordinary one- 
dimensional fermions in the quantum-mechanical analogy (with X 2 playing the role 
of imaginary time). An identical joint distribution arises also in the study of one- 
dimensional non-intersecting random walkers |19|. 

Returning to the case of fixed V, and still considering r ~ N, we note that the 
N 2 — N "gauge degrees of freedom" of V enter only through N combinations V aa . An 
application of the central limit theorem then implies that, after integrating the "gauge 
degrees of freedom" out, V aa become independently distributed with the Gaussian 
weight exp {— N 2 V 2 a /2TiV 2 ^ . Consequently, the form of the kernel coincides with 
the case when V is random. 

Let us turn now to the case when r ~ 0(1). The remaining step in the derivation 
of the joint distribution function of the eigenvalues is to implement the integration 
over dfir r -\ ({a re }). This can be done in two ways which lead to different forms of 
the distribution functions but are ultimately equivalent in the N — > oo limit when 
applied to the calculation of the correlation functions. The first - exact - procedure 
requires enforcing the conditions of orthonormality between the r complex vectors a K . 
Although quite simple in the r = 1 case, already for r = 2 it leads to a somewhat 
lengthy calculation, as we will show below, and it quickly becomes unmanageable for 
larger values of r. If, however, the ultimate interest is in the correlation functions R nm 
in which the scaling limit N — > oo is taken at fixed values of n and m, integrating out 
the remaining 2N — n — m "energies" leaves the correlation function dependent only 
on a finite number of components a Ka of the r vectors a K . According to a well-known 
result in RMT 0, the latter are distributed independently according to the Porter- 
Thomas formula e _JV ' a " Q l . The distribution functions obtained by replacing the whole 
of dfiM ({a K }) with the product of independent Porter-Thomas distributions for each 
of the Nr components of a K can be interpreted as resulting from the substitution of a 
"canonical" ensemble for the gas of levels in place of the "micro-canonical" ensemble 
implied by the identity Trif' = Tr(H + V). 

2.1. Exact ("micro-canonical") joint distribution functions 

To illustrate the method, let us first consider the r = 1 case. The corresponding 
distribution function and the lowest-order nontrivial correlation function (i?n in our 
notation) have been derived using a somewhat different approach in Ref. V is 

parametrised as vNsl ® a^ , and for a single normalised complex vector a the measure 
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dfxn) (a) becomes 

dA*(i) (a) = JiS (a f a - l) e-" ata d [a] (17) 



where d [a] denotes the unconstrained measure n^Li da a da*/7r, and 77 is a positive 
infinitesimal introduced for convergence. The normalisation constant J\ is determined 
from 

1 = Ji /d W r ge*-*-*^- = Ji r ^-^-iv = T^W" (18) 
J J -00 2tt J-oo 27r (ia: + 77) (A - 1)! 

Substituting this measure into l|16|l we obtain 

C v \c\ 2 fi N ^_ t n iW(ea)A . (2V-1)! 



^ ({<}, {*'}) = ^T^e- E2 - W(£a) A (A) A (H 



(2nf 2 - N " ' " '(N\v\)»-i 

XS I J2 e '{3 ~ J2 €a ~ NV } det0 [ Signl1 ( £ /5 ~ £q )] ' ( 19 ) 

As expected [22 1 the determinant is nonzero only if the two sets of levels satisfy the 
"interleaved combs" constraint e a < e' a < e Q +i for v > and e a -i < e' a < e a for 
v < 0. 

We next consider the r = 2 case. It is convenient to use a symmetrised expression 
V = vN (a (£> at — b ® bt) where a and b are two complex mutually orthogonal N- 
component vectors of length 1 , and v is the scalar (positive) parameter characterising 
the magnitude of V. More general expressions can be reduced to this form by a simple 
rescaling in the final results. It is also convenient to re-exponentiate the (5-functions 
in the determinant in (|16fl so that the A integrals are performed at a later stage in the 
calculation. The measure d/i( 2 ) ( a , b) is explicitly parametrised as 

dn {2) (a, b) = J 2 d [a] d [b] 6 (a, b) e -v(^*+ h ' h ) (20) 
where 8 (a, b) stands for 



S (at a - 1) S (b*b -1)8 (a*b + b*a) 5 Q [a+b - b+a] ^ , 



(21) 



and Ji is the normalisation constant. To compute J 2 , the delta-functions are 
exponentiated with the help of two real Lagrange multipliers x and y, and a complex 
Lagrange multiplier z: 

1 = J 2 / [da] [db] (5 (a, b) exp {-77 (a*a + b*b) } = J 2 / dxd y Mz ^ z [ da ] [db] 
J J (27r) 

x exp {i (a; — iry) (l — a' a) + i(y — i-q) (l — b^b) — iza^b— iz*\J&} . 

Performing the integrals, we find 



2N 



J 2 = 47r^-J (-1)^ (jV - 1)! (AT - 2)! (22) 
With these definitions, we have 

r i r N 

V ({e} , {e'» = — ~~]fjrA (H) A (#') / d// (U x ) JJ dA fc d [a] d [b] 5 (a, b) exp {*Tr (XV)} 
(2tt) 7 fc=1 



x (- 1 )' T+ ' T ex P S * (eaA CT(Q ) - e^A CT / (Q) ) \ exp ^ - ^ W (e a ) ^ 



(23) 
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Proceeding as in the computation of J2, we obtain 



V({e},W}) = 



Cu I J2I |c| UN 



(27T) 



N 2 



A(H)A(H>) (jLy J dxdydMzd^z Sx+y) 



N 



fc=i 



exp 



|iA fc (e 



/Nv\ 



N 



(A fc - x + iry) (A* + y - 477) + |z|' 



a=l 



After integrating out the eigenvalues of A we obtain 



T ({«}.{«'}) = 



Cu|J 2 ||c|VjvJV! 



(2tt) 



N 2 -N 



A (ff) A (#') 



2^q=1 \ e " e cr(a 



dr exp • 



2Nv 



(t 2 -C 2 )(-^) 



A' 



t 2 -C 2 



iV 




+ 9 (C 2 - i 2 ) 



where t = x — y, t = x + y, £ = 2\z\, s t = Signi, the ambiguities in the 
analytic properties of the integrand are resolved by letting t — > t — irj, and A a 

2^Q = 1 



c <r(a) 



/2Nv. The extra factor AH comes from performing one of the sums 

over permutations. Note that the parameter A a depends only on the absolute values 
of the distances between the levels in contrast to the case of rank r = 1 perturbation 

EP 

Integration over r gives a factor of AttNvS \Y) a (e a ~ e ' a )] which simply enforces 
the identity Tr H = Tr H' . Performing the remaining integrations, we find 



2Cu \c\ (in 



(27T) 



N2 - N \Nv 



N-l 



— N\(N-1)\S 



xA(H)A(H')e 



£(-i) ct 0(i-a ct )(i-a ct ; 



N-2 



(25) 



Note the hard constraint on the sum of the absolute values of the level shifts 



y 



V(a) 



< 2Nv. Note also that this result cannot be simplified in an obvious 



way in the N — > 00 limit since A CT does not scale as l/N. Indeed, the natural scale for 



v is the mean level spacing A, and ^ Q 



<V(a) 



O(NA) so that the A CT - 0(1). 



The distribution functions (|16f) and (|27|l below are equivalent only in the sense that 
they lead to identical correlation functions in the "thermodynamic" limit N — > 00. 

2.2. "Canonical" joint distributions 

As discussed above, after neglecting 0(1/N) corrections, the measure d/Z( r ) for finite 
r is brought to the simple form 

r N 

l e -"l*-l". (26) 
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Due to the factorised structure of this measure, the transition kernel in acquires 
the determinantal form 

detV7 £ ] S(e a -e p ) (27) 

ap Lai j 

where the action of the differential operator 2?{ c; i} on an arbitrary function /(e) is 
defined as 

T> {e - M f(e)=det(i-M/de)f(e) (28) 
and its inverse has a convenient integral representation 

/n^^exp{-M 2 }/( e + ^). (29) 

K=l 7F 

In the opposite case of infinite rank, the central limit theorem ensures that 
the diagonal matrix elements V aa are independently Gaussian distributed. In 
consequence, the differential operator £>{ e; {)} simplifies to 

P {£;C} =exp|-i T rf) 2 ^|. (30) 

Here and below when considering the infinite rank case we assume Tr V = since a 
non-zero value of Tr V can be accounted for by a trivial uniform shift of the spectrum. 

3. Parametric correlation functions 

3.1. Parametric correlation functions in the energy representation 

The objective of this section is to present a derivation of the multipoint correlation 
function R nm [equations (|38|1 and (|44|) below]. The joint distribution function defined 
by equations Ijltjfl and l|27|) serves as the basis for this calculation which employs 
a scheme based on the orthogonal polynomial expansion. In accordance with the 
definition the parametric correlation function is represented as 

N N 



(31) 



Rnm ({EaYU » = / f[ de a f[ de' V ({e} , {e'» S({E})S'({E'}), 

J a=l p=l 

where the source terms S({E}) and S'({E'}) are defined as 

71 N m N 

™-nE^ E «- e «.)' s '({ E '}) = T[Y, 6 ( E ' b -tk)- (32) 

a=la a =l b=l p b =l 

Integrating over all e' we obtain 



({E a } n a=1 ,{E' b }Zi)= [ f[de a A(e)exp\-Y,W(e a )\s({E}) 

a=l { ct=\ ) 



N 

x TT vj 1 n 

(3=1 



m N 

A(e)nE^- e ^) • ( 33 ) 

6=1 /3 b = l 

Let us denote as 7r n , and 7r m , respectively, the set of all distinct indices a a in the 
multiple sum over a±, . . . ,a n , and the set of all distinct indices (5 b in the multiple sum 
over /3 m . It is now convenient to split the product of 2? _1 's as 

N 

ii *r«u n **u- (34) 
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All the factors in the second product act only on the variables in A(e) while the 
operators in the first product act on the variables in both A and the source term 
delta-functions. The Vandermonde determinant involves only polynomials of degree 
p— 1 in its p-th row, so that the action of a product of N identical differential operators 
transforms the elements in its p-th row into different polynomials of the same degree 
and with the same leading term. Consequently, 

JV N 

Uv^A^^Aie), and JJ T>^ } A(e) = JJ %,;,}A(e). (35) 
Integrating by parts over the variables in w m we arrive at 



6=1/ 



11 {<W> 




6=1 



J2w(e a )\A(e)Y[5(E a -e c 



(36) 



where t> 



r^-. — is{ e ._{,}. The above expression is manifestly asymmetric with respect 
to the confining potential W(e). As long as all the energies are close to the centre of 
the band, and the change in the average density of states induced by the potential 
V is negligible, the action of T> on W can be ignored. It is instructive, however, to 
consider a more general case. To this end, we introduce a new operator 



V {t . M = exp {-W (e) /2} V {t , i} exp {W (e) /2} . 



(37) 



When H36fl is rewritten in terms of T> and T> , each of the Vandermonde determinants 
appears in combination with I1q=i ex P {—W (eq) /2}, so that they can be equivalently 
written as det Q /3 ip a -i(e/3) where ip a -i(e) = exp {— W (e) /2}p Q _i(e), and p tt (e) are 
the polynomials orthogonal with respect to the weight exp{— W(e)}. Expanding the 
determinants, and using the orthonormality of the set of functions (p a (e), we arrive 
at the following representation of the correlation function as an (n + m) x (n + m) 
determinant: 



({Ea} n a=1 ,{E' b }Zl 
( 



det 



Wei 



V( E , [k(E ax ,E' b2 )-S(E ax -E' b2 )] 



V^V^KE'^El) 
where the indices a and a 1 (b and b') run from 1 to n (from 1 to m). The kernel 



(38) 



N-l 

k(E,E') = Y, <Pv{E)^ P {E') 

p=0 



(39) 



describes the correlations among the levels of the original matrices H, while the lower 
right corner of the matrix describes the correlation functions within the ensemble 
defined by the distribution function V (H) = Ctj exp {— Tr W (H — V)}. The cross- 
correlations between H and H 1 = H + V are encoded in the off-diagonal blocks of 
(j2HJ- In the bulk scaling limit the kernel k takes the form JT] 

k(E,E') = 3 exp {inujp (e)} /u = np (e) k [irp (e) u] , (40) 
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where oj = E' - E and e = (£" + E) /2. 

The bulk of the present paper is devoted to the case when the strength of the 
perturbation V is such that the parametric correlators between the "old" and the 
"new" levels are 0(1). This condition is fulfilled if TrF 2 ~ 0{N 2 A 2 ). I n this regime 

the operator V^.^y can be explicitly evaluated using a representation of T> as an 
integral over a vector of Grassmann variables X : 

V { ,. M ee exp [-W (e) /2] J d X d X * exp j- \ X \ 2 - X*v X j^ exp [W (e) /2] 

~ J d X d X * exp [W' (e) X *v X /2] exp j- \ X \ 2 - X*VX^ 

= det[l-vW'(e)/2]T> {e;im , (41) 

where i(v) = v/ (1 — W (e) /2). Assuming that all the energies in both sets {E} 
and {£"} lie within a few level spacings of some central energy e, and since the 
determinantal structure of Ij38(l ensures that the operators V and X> _1 come in pairs, 
the determinant in the prefactor can be dropped. 

The derivative of the confining potential is related to the real part of the average 
Green function via 

W (e) = 2l ee 2m (e) . (42) 

Consequently r (v) can be represented as the solution of the equation 

r (v) = v + v?StG (e) r (v) , (43) 

which can be immediately identified as the equation for the reactance matrix for 
scattering off the potential V. In fact, as can be shown using the methods of statistical 
field theory |12j . the parametrisation of T> in terms of the corresponding reactance 
matrix retains its validity in the more general case of non-invariant distributions 
V(H). The eigenvalues of the reactance matrix are conveniently parametrised as 
tani5 K where 5 K are the corresponding phase shifts. 

As discussed above, the parametric correlations generically fall into two regimes 
corresponding to the rank of V being finite or infinite. In the latter case, the condition 
TrV 2 ~ 0(7V 2 A 2 ) ensures that, for a generic eigenvalue Vi of v, and at a generic 

position e in the spectrum, ViW'(e) <C 1. Thus the distinction between V and V is 
meaningful only in the finite rank case, and, to simplify notations, we will use the 
symbol T> in both cases with its parametrisation implied by the context. 

Similarly, to the leading order in the 1/N expansion, the density of states p (e) can 
be approximated by an e- independent inverse mean level spacing 1/A, thus simplifying 
ff to the form quoted in |17| : 

%im {{Ea} n a= l ) { E 'b\ 



b=l) 



k(e a2 -e ai ) Vr l ,(k[e' b -e ai ] ~S[e' b -e ai ]) 
det I K ;v ) I , (44) 

V {e a2 ;v} k ( £ a 2 -S' bl ) fc^-^J 

where v = v/A, and U\ nm — A n+rn R nm . In the r — > 00 limit this equation reproduces 
the results obtained previously Ej ■ The representation of R nm hi the form of A 
integrals characteristic of the non-linear sigma model calculations ^1 corresponds 
to the Fourier decomposition of Vk and V^ 1 [k — 8] in Q44fl. 
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3.2. Universality in shifted distributions 

To draw a connection with some previous studies of parametrically deformed random 
matrix ensembles ESI , let us briefly consider the opposite case of a perturbation 
strong enough to significantly affect the density of states. In this regime, the cross- 
correlations are small as 1/N. There remains, however, the question of whether 
the correlation functions of a matrix ensemble defined by the shifted distribution 
Vv (H) = Cjj cxp {— Tr W (H — V)} are universal. To the best of our knowledge, 
this question has been so far answered in the affirmative only in the case when W is 
quadratic (23 ( see a ls° HI] where the case of non-Gaussian potential with a linear 
source term was considered). The formalism developed above affords an opportunity 
to extend the proof to the case of arbitrary W(H) (still, however, assuming strong 
confinement). The corresponding correlation functions, as follows from (|38f) . have 

the determinantal structure with the kernel k„(w,e) = T>^ ti .~^D^ e2 .^\ (w, e). The 

operator T> in this regime is given by equation (|3L)fl . Expanding the difference 
W (e 2 ) - W (ei) m luW (e), we obtain 

Ww,e) = exp{-LuW (e) /2} exp {-Tr 2 x 2 d UJ dA exp{uW' (e) /2} ^ exp ^ lu}7r P ( e )> (45) 

TTU) 

1 11 . . ... 

where x = (TrO 2 ) /it. By virtue of l|42(l . the last two exponents can be combined to 
read 3 exp {ojQ (e)}, where (5 (z) is defined as the analytic continuation of the retarded 
Green function G(e) = j- dzp (z) / (z — e) + iirp (e) onto the complex plane with cuts 
along the real axis outside the support of the density of states p (e) . By employing an 
integral representation of the differential operator 

exp{-Tr 2 x 2 d u d e } = J ^^exp{-|C| 2 A 2 ^-C*^ + C9 e },(46) 

and differentiating with respect to e, equation (|45|) is brought to the form 
d e [e^{coW'{e)/2}k v (co,e)] 

{-ICI 2 /ttV} 3exp{(c - C) Q (e + Q}G' (e + 0.7) 

After a change of variables r — \(\ 2 , u — e larg< ' we note that, since Q is analytic 
by construction, the only non-analyticity of the integral over u is due to the pole 
generated by the solution u (e) of the equation n 2 x 2 G (e + u) + u = (existence of 
multiple solutions of this equation would contradict level number conservation). We 
thus find, upon rescaling u and integrating over r, 

d ( [exp {ojW (e) /2) k v (cu, e)] = 9f / du ^^^^L g' (e + u) 

(48) 



d(d(* 
— T-T ex P 



tt 2 x 2 G (e + ?i) + u' 
exp {-wit (e) /7r 2 a; 2 } 5' (e + u (e)) 



1 + ttW (c + u (c)) 
where the contribution from integrals along the cuts vanishes upon tak- 
ing the imaginary part of this expression. Integrating over e and using 
[1 + TT 2 x 2 g' (e + u (e))] (rfuo/rfe) - -^ 2 x 2 ^' (e + u (e)), we find 

kJw, e) = exp {-^^ (e) /2 + triShto (e)} (£) , (49) 

w 

where the renormalised density of states is p v (e) = — *<su (e) /tt 2 x 2 . The exponential 
prefactors cancel out due to the determinantal structure of the correlation functions, 
thus establishing their universal form. 
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3. 3. Parametric correlation functions in the number representation 

The "complementary" number-dependent correlation functions are generically 
expressed as derivatives of the corresponding generating functions. The archetypal 
generating function is the probability ^ nr! /(J, J') to find n(J) = n levels in the 
interval J of the unperturbed sequence and n' levels in the interval J' of the perturbed 
sequence, where n( J) is the number of levels in the interval J: 

^ n ,(J,j' ;? ) = ( W) 5 n , >n(JO > = ^ I - 7r ^ E (fc r n )!(fc'-n')! - (50) 

Here r^fc' represents the correlation function Rkk 1 integrated over the direct product of 
intervals J k ®J lk with the corresponding measures dmj and dmji . The intervals J and 
J' may, generically, consist of an arbitrary number of disjointed segments. Equation 
H5Ufl is a straightforward generalisation of the non-parametric analogue P n {J) (see, 
e. g., Appendix 7 of Ref. [2]). For later convenience in carrying out the summations, 
r^k' can be represented in the form of a fermionic functional integral 

r kk , = detK [ V^V^ ( f dmj(u)Z(u)£(u)) ( f dmj>(w)fj(w)ri(w) 



xexpjy du J dwV(u)K 1 (u,w)i& (w)j (51) 

where ^ = (£, fj) is a fermionic doublet. The matrix kernel K is easily reconstructed 
from (PH>: 

K ^={^tr\ (52) 

V Vlk{u,w) V w L Vlk(u,w) J 

the matrix elements of the scalar kernel k are given by l|39(l . and I is adopted as a 
shorthand for the Dirac 5- function. For brevity we also drop the subscript x(v), and 
V is understood to be the function of cither the single parameter Tr v 2 or the r phase 
shifts 8 K parameterising r (v) in the infinite, and finite rank cases, respectively. 

Equations l|5Ufl . (|51|l . and (|52|) can be used directly to compute a certain class 
of correlation functions exemplified, e. g., by the correlation between level spacings 
of the original and perturbed level sequences centred at given energies E, and E' , 
respectively. Since our main interest, however, is in the correlation functions in the 
number representation, we need to be able to relate the levels in the new sequence 
to the parametric "descendant" of a given level e a in the unperturbed sequence. The 
absence of level crossings allows to do this by demanding, e.g., that if an interval 
J e = (— oo, e] of the unperturbed sequence contains exactly n levels (as counted from 
the bottom of the band), the corresponding interval J^ +Ul of the perturbed sequence 
contains exactly n + q levels, where q is fixed and n is arbitrary. By construction, 
then, the uppermost level in the interval J' e+U is separated by q — 1 levels from the 
parametric descendant of the uppermost level in the interval J e . Summing over all n 
we find the corresponding probability as 

00 ^ 00 / ^k+k'+q ^ 

y q (J e ,J' e+u ;x) = (*„ (J .),»(j^ M ) +fl ) = E n \( n + q )\ E (*_„)!(*'_„*%)! 

n— k—n 

k' —n-\-q 

d± i^ y rfcfc ' 4^', (53) 

2tt ^ k\k'\ + y ' 

k=0 
k' —q 
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where z± = 1 + e ±l ^. Substituting l|51|l into i|53|) . one finds 



% q ( J, J'; r) = (-l) ? ^ 27r (w, y; 



where 

Z 9 <£) = det 



g— 1 ^fc' 

lao-Xn^)] - £ -^det 



fc'=0 



Icr - KILy ((/)) 



(54) 



(55) 



7=0 



and, denoting the Pauli matrices as at, 



n(0) = 



o 



o 



n 7 (0) 



o 



o 



The determinants are understood as functional determinants on the space of two- 
component functions defined on the product interval J ® J'. 

Although the formalism in principle can be applied at the spectral edge, we are 
primarily interested in the correlation functions in the bulk of the spectrum. Taking 
the scaling limit would then greatly facilitate the calculation. In order to avoid 
unessential complications related to the energy dependence of the average density 
of states, we replace k with its value (|40|l in the scaling limit assuming A = p _1 (e) is 
a constant, and then regularise k as 



k n (u — w) = 



_ sin 7T (m- w) (1/2MM+M) 



(56) 



7T (u — w) 

where the limit jy — >■ is implied in all expressions involving this kernel, and the 
energy variables are rescaled by A. Using equation l|51|) it is easily shown that 
(n(Joo)) = <n(J^)> = 2/r?, and ([n(J M ) - n(J^)] 2 ) ~ 0( V ). Thus, although the 
regularisation formally violates the level number conservation, the corresponding error 
tends to zero in the limit 77 — > 0. At the same time, for any finite value of rj the number 
of levels in the semi-infinite intervals J e and J' e+W is finite, ensuring the convergence 
of the integrals over energies at the lower limit. In the following we will suppress the 
index 77. 

As written, l|54|) involves a matrix oscillating integral kernel defined on a product 
of semi-infinite intervals. However, as we will now show for the case q = 0, it can 
be rewritten in the form which is (i) more amenable to numerical analysis, and (ii) 
makes the integrability of the kernel (in the sense discussed in Ref. [2HJ) manifest. 
Without loss of generality, we can set e = 0, and shift the variables so as to define the 
determinant on the quadrant (— 00, 0] <E> (—00, 0]. The corresponding shift operator is 
absorbed into the redefinition V — » e ^ d / du f). The term involving the S- function in 
the upper right corner of l|52")) can be separated to reveal the dyadic structure of the 
remainder: 



K 



V 



( k T>- X k 



V- 1 ! 




where the we have also used 2?J u fc(ii — w) = t> u k(u — w). Now, using the identities 
(Icr + z-T>- 1 1<t + )- 1 = (loo - z-V- 1 Ia + ), and det(Icr + Z-T)- X la + ) = 1, where 
(7+ = (o"i + z<72)/2, one obtains 



det 



a Q -KU((t)) =dct I- (z + k - z+z-k (v^lj T> + z- (v^k^Vj .(57) 
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Employing the Fourier representations, 



fc(tt) 
S(u) 



dX 



A)0(1 + A) 
1 



one finds 



z + k - z+z-k [V^ 1 !] V + z_ [V^k) V (u,w) 



2m 



dix—^ + z-V- 1 (A) V \, (58) 



A 



iS 



where in the infinite rank case 

V (A) = exp {inuX + tt 2 x 2 X 2 /2} , 
and in the finite rank case 



V (A) = exp {iiruj\} Y[ [1 — iA tan Sk] 



(59) 
(60) 



Finally, the cyclic invariance of the determinant and the identity z + + z_ = arc 
used to perform the integrals in the u-w space, bringing the determinant to the form 



Z q {<jJ,tA) =det lcr - ATI (</>) = dct [I — /C (</>)] , 



where the matrix elements of the kernel K, (<fi) are 



K, (A, /i; 



= — V^(a)p(m) 



F(A;0)-FO*;0) 



and 



F(A;< 



_ (g++^) r d ^M 

iri J-oo fi-X 



(z + -z-)V-\X) 



(61) 



(62) 



(63) 



Here the integral is understood in the sense of the Cauchy principal value, with 
variables A and \i restricted to the interval [—1,1]. We employ the notation y to 
denote the parametric dependence of Z in both finite and infinite rank cases. A 
similar, although somewhat more cumbersome, expression, is easily reconstructed for 
the contribution of det \oq — A±I 7 (<fi) in the q ^ case. 

Equations l|61(l . 1)62(1 . and 1|63|) form the central results of this section. In the 
remainder we will consider several applications of these equations to the calculation of 
number-dependent correlation functions. The integral of the generating function itself 
J dcfiZo (lu, y; </>) /27r gives the probability ^Po(w, y) = ^Po (Je, Je+ui)- In order to obtain 
the correlation functions P q and C q introduced above, we rewrite their definitions in 
the form 

P g M = (d e n(J e )d e n(J^)S n{jMj , + ^ )+q ), (64) 

C q (?) = (d r n (J e ) d r n (J' e+U ) S nijMj , +a)+q ) . (65) 

Using translational invariance both in the energy and the parameter spaces, rescaling 
the energies by A, and rescaling the variables in y by C 0, we obtain for the 
rescaled correlation functions p q and c q j2H] 



2n 



2tt z + Z- 

,-2-k 

dw 



(-dZ)z q (u,r, 



2tt z + Z- 



(-5 r ) Z 9 (w,y; 



(66) 
(67) 



Parametric Spectral Statistics 



16 



Let us now consider several special cases, (i) In the infinite rank case y reduces to a 

1 /2 

single parameter x = (Trv 2 ) /tt. The single level universal velocity autocorrelation 
function c (x) = cq (x) is plotted in Fig. ^ For comparison we also show the 
results of direct numerical simulation of large random matrices |27| . We also find, 
based on numerical evaluation of equations (|f>f>|t . 1)61(1 . ((62(1 . an d 1)63(1 that the 
function p (to, x) = po (w, x) - the distribution of parametric single level shifts - has a 
Gaussian form p (u), x) = exp in accordance with the earlier 

conjecture [2D1 - An analytical proof of this statement, as well as finding an analytical 
expression for the width a (x) (shown in the inset in Fig. ^) are some of the open 
questions posed by the results of the present investigation. 
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Figure 1. Single level universal velocity autocorrelation function obtained from 
1671 . 1611 . 1621 . and 1631 (solid line) vs. direct numerical simulation of large 
random matrices [2TJ (dots). The width a (x) of the Gaussian distribution of 
level shifts together with the x 2 asymptotics at small x is shown in the inset. 



(ii) In the finite rank case further progress can be made due to a simple pole 
structure of V^ 1 (A). Leaving a detailed consideration of the finite rank case to a 
future publication, we present here the results for Co and po in the r = 1 case, y 
is now parametrised by a single phase shift 5. In Ref. ^2] the distribution po for 
r = 1 was found using the fact that, as a special property of r = 1, this distribution 
coincides with the level spacing of the combined ensemble of levels {e a } U {e^}. Using 
the present formalism, this quantity, as well as the corresponding cq can be evaluated 
in a more direct way. Introducing an auxiliary non-parametric kernel 

6 ( A)M ) = -L 1 - exp f^ A -^ i (68) 
its Green function g (A, /i) = (I — 6) , and the determinant d(ui) — det (I — I), we 



find (assuming S > 0) 

Cn r_1 ' (5) = — -pr f — d (w) exp {— ttw arctan<5} / 
dd z J 2vr 1 J_ 



arctan o — iX 
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A somewhat more lengthy expression is obtained for p : 

9 (\ /-0 exp {iirujn} 



(r=i) f X \ or s^d(uj) f 1 g (A, 

Po ' [oj; o) = a (00) — - — exp { — irui arctano} / dXd/j, - 

2 J -i 1 

1(1 — i/itan5) 



— iXt&nS 

+ (1 -i/ztan<5)i(w) -T{uj,5) \ , (70) 



tan S 

where t (a>) = J_, dXdftg (A, /i) exp {i7ro; (/i — A)}, and 

T(oj,<5)= / dXd/jg (A, /x) exp {i7ro; (fi — A)} (1 — i/itan5) . (71) 



Up to a simple transformation, the kernel t is identical to the kernel which controls 
the non-parametric level spacing distribution ©x- I n light of the "interleaved combs" 
constraint on the evolution of levels when r — 1, such a connection to Si is not 
unexpected. 

Within the framework presented here various number-dependent correlation 
functions in the unitary ensembles can be computed "at will". Nevertheless, several 
open questions remain unanswered, and their resolution would greatly enhance our 
understanding of parametric evolution of random matrices. As already mentioned, the 
Gaussian form of po in the r — ► oo case is probably an indicator that the <j> integration 
can be performed analytically in some form in the infinite rank case as well. If so, 
a more "processed" expression for Cq^°° may emerge, throwing more light on the 
properties of the kernel K.{4>). Separately, an analysis of the asymptotic behaviour 
of the generating functions along the lines of Ref. [2Hj is a natural way to develop 
the present formalism further. Finally, the technique of "counting" the number of 
levels between the spectral edge and a point in the bulk of the spectrum by means 
of introducing a regularised kernel should be applicable to ensembles with different 
symmetries, so that if an analogue of l|38|l is obtained for such ensembles, it is likely to 
open the door to the computation of the corresponding number-dependent correlation 
functions. 
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